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Abstract:  We  discuss  the  problem  of  computing  the  intersection  curve  of 
two  algebraic  surfaces,  each  of  which  possesses  rational 
parametrization.  The  special  case  where  the  two  surfaces  are  quadric 
is  analyzed  in  detail,  using  a  general  decomposition  theorem  which 
guarantees  the  existence  of  a  simultaneous  canonical  reduction  of  two 
quadratic  forms  in  Euclidean  n-space.  In  homogeneous  4-space  this 
yields  a  classification  and  simple  parametrization  of  all  possible 
intersections  between  two  quadric  surfaces.  Using  these  results,  we 
treat  the  problem  of  analyzing  the  structure  of  solid  bodies  defined  by 
Boolean  combinations  of  half-spaces  bounded  by  arbitrary  quadric 
surfaces.  The  analysis  given  leads  to  fast  versions  of  some  of  the 
procedures  required  for  geometric  modeling  systems  which  admit  general 
quadric  surfaces  into  their  vocabulary  of  basic  shapes. 
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1.  Introduction. 

In  his  survey  [We80]  of  geometric  modeling  systems,  M.  Wesley 
makes  the  following  cogent  remarks:  'There  are  some  fundamental 
problems  to  be  solved  in  the  field  of  geometric  representation. 
Firstly,  curved  surfaces  are  required;  in  the  short  term,  cylinders  are 
probably  sufficient  for  many  applications,  to  be  followed  by  cones  and 
then  more  general  curved  surfaces.  Secondly,  the  geometric 
representation  has  to  be  absolutely  robust  under  all  conditions.  The 
system  can  be  expected  to  be  used  to  generate  arbitrary  objects  without 
human  intervention.  It  is  essential  that  the  objects  always  be 
numerically  and  topologically  valid;  e.g.  it  must  handle  situations 
where  surface  details  become  arbitrarily  small  without  generating 
invalid  topologies.  Thirdly,  the  system  must  allow  computation  of 
necessary  results  with  optimum  efficiency  in  terms  of  computation  time 
and  storage  efficiency.  As  the  applications  leave  the  stylised  blocks 
world  of  the  laboratory  and  approach  the  full  complexity  of  real 
engineering  situations,  the  limits  of  computational  capacity  can  be 
expected  to  be  reached,  and  choices  between  modeling  systems  may  be 
based  on  cost  rather  than  the  ability  to  perform  the  operation  at  all.' 

If  for  the  moment  we  ignore  the  computational  cost  factor  on  which 
Wesley  lays  a  very  justified  stress,  it  becomes  possible  for  an  ideal 
geometric  modeling  system  to  deal  with  arbitrary  bodies  and  surfaces. 
There  do  exist  general  algorithmic  techniques  for  deciding  all  the 
questions  about  such  bodies  which  a  geometric  modeling  system  needs  to 
ask.  These  techniques  derive  from  Tarski's  general  method  for  deciding 
quantified  formulae  in  the  elementary  theory  of  reals,  which  has  been 
improved  substantially  by  various  other  authors  in  the  years  since 
Tarski's  original  paper:  see  [Ta51 ,Co75,Ar81 ,La77 ]  and  [SS81]  for  a 
review  of  some  of  this  material. 

However,  these  general  methods,  which  are  based  upon  systematic 
but  unoptimi^.ed  use  of  elimination  theory  and  of  Sturm's  location 
theorem  for  the  roots  of  real  polynomials,  are  too  expensive 
computationally   for  use  in  a  practical  geometric  modeling  system.   For 
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thls  reason,  one  must  focus  on  a  class  of  surfaces  whose  Boolean 
combinations  can  be  handled  with  considerably  greater  efficiency  than 
can  algebraic  surfaces  in  general.  The  following  definition  specifies 
a  class  of  surfaces  which  is  advantageous  (but  not  necessarily 
sufficiently  advantageous)  in  this  regard. 

Definition  1:  A  surface  S  is  said  to  be  rational  if 

(a)  S  is  defined  by  a  polynomial  equation  P(x,y,z)  =  0  having  rational 
(or  at  least  real  algebraic)  coefficients; 

(b)  There  exists  a  parametrization  x  =  X(u,v),  y  =  Y(u,v),  z  =  Z(u,v) 
of  the  surface  S  by  functions  X,Y,Z  which  are  quotients  of  polynomials 
in  u,v  having  rational  (or  at  least  real  algebraic)  coefficients.  (We 
suppose  this  parametrization  to  be  defined  for  u,v  in  some  domain  D 
which  is  itself  defined  in  a  simple  way  by  polynomial  inequalities  in 
u,v,  and  suppose  that  except  on  a  few  equally  simple  curves  and  points 
this  parametrization  is  1-1.) 

The  advantage  of  working  with  such  surfaces  is  that  their 
intersections  can  be  deterrained_  relatively  efficiently,  since  if  Pj^  =  0 
is  the  polynomial  equation  for  one  such  surface  S,  and  X9,Y9,Z9  give 
the  parametric  representation  for  another  surface  S-,,  then  the 
intersection  Sj^*So  is  represented  by  the  set  of  pairs  satisfying  the 
simple  equation  P  j^  (X2(u,  v)  ,Y2(u,  v)  ,Z2(u,  v) )  =  0.  Although  a  related 
algebraic  condition  can  be  written  even  if  the  surfaces  S  ]^  and  S2  are 
not  rational  in  the  sense  defined  above,  resultants  would  have  to  be 
used  to  derive  this  condition,  substantially  degrading  the  efficiency 
of  the  calculation  required  to  deal  with  S]^*S2. 

The  class  of  real  rational  surfaces  is  extensive.  Examples  are  as 
follows.    Any   surface  of  the  form  z  =  P(x,y),  where  P  is  a  polynomial 

or  a  quotient  of  polynomials  is  evidently  rational.   This  includes   the 

999 
paraboloid   and   the  hyperboloid  of  one  sheet.   The  sphere  x'"+y"'+z   =  1 

is  rational,  since  it  has  the  (stereographic)  parametrization 
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f  N    /  1-u  -v"     2u       2v   . 

(x,y,z)  =  (   ^   ^.  — T-^.  — ^^--f) 
1+u^+v-   l-Ki'-H-v"   1+U-+V" 

9   9   T 

Similarly,  the  cone   ^--x^-y"   =   0   is   rational,   since   it   has   the 
parame t r i z a t i on 

(x,y,z)  =  v( _,  _,  1) 

1+u-   l+u- 

9    9 

The  cylinder  x^+y   =  1  is  rational  and  has  the  parametrization 

o 

(x,y,z)  =  ( _,  _,  v) 

1+u-   l+u- 

999 
The  double-sheeted  hyperboloid  z  -x^-y''  =  1  has  the  parametrization 

9   9 
,      ^    ,   2u       2v     l+u^+v". 
(x.y.z)  =  (     ,.    ,  ,,    ,  ,) 

l-u^'-v'"   l-u^-V"   1-U--V"' 

where   u   and   v   vary   in   the   open  unit  circle.   The  single-sheeted 
hyperboloid  x"+y  -z~  =  1  is  a  ruled  surface  and  has  the  parametrization 

,      .    ,l-u^+2ut   2u-t(l-u'^) 
(x,y,z)  =  ( — ,  ^— i,  t). 

1+u-       1+u^ 

The  torus  formed  by  sweeping  a  circle   of   radius   r   <  R   and   center 

909 
(R,0,0)  around  the  circle  x~+y-=R'-  in  the  x-y  plane  has  the  equation 


or 


|(x,y,z)   -  R_liilZli^|2   =   ^2 
(x2+y2)l/2' 

(x2+y2+z2+R2_^2)2   =   4R(x2+y2), 


and  the  parametrization 

(x,y,z)  =  (ri:4^R).(il4.  A->  0)  +  (^>  0,  -^). 
1+v-      1+u-   1+u-  l+v2 

Many  other  surfaces  are  rational  also.  For  example,  any  surface  of 
revolution  of  a  pararaetrizable  curve  (x,y,z)  =  (f (t ) ,0,g(t ) )  about  the 
z-axis  is  rational,   and  has  the  parametrization 

(x.y.z)  =  (f(t)ll^.  f(t)Jl-,  g(t)). 
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and  a  defining  equation  which  can  easily  be  obtained  by  elimination 
from  the  parametric  representation  of  the  initial  curve.  Note  also 
that  any  affine  or  projective  image  of  a  rational  surface  is  also 
rational. 


2.   Intersection   of   Quadric   Surfaces ;   Algebraic   Analysis   of    the 
n-dimensional  Case 

It  is  easy  to  see  from  the  above  remarks  that  all  real  quadric 
surfaces  are  rational.,  Indeed,  if  we  proceed  protectively,  writing  the 
polynomial  equations  for  all  our  surfaces  as  homogeneous  polynomial 
equations  in  four  variables,  then  the  equation  for  any  quaciric  can  be 
written  as  La« a  =  0,  where  L  is  a  4x4  symmetric  matrix,  and  where  a  is 
a  4-vector  (x,y,z,w).  We  can  diagonalize  L,  which  corresponds  to 
carrying  out  an  appropriate  projective  transformation  of  our  surface, 
and  normalize  the  transformed  L  so  that  each  of  its  diagonal  entries  is 
either  ±1  or  0.  Then  this  normalized  L  can  be  brought  to  a  form  in 
which  it  has  at  least  as  many  positive  as  negative  entries,  and  in 
which  the  diagonal  entries  of  L  are  arranged  so  that  the  ;^ero  entries 
precede  the  positive  ones  which  in  turn  precede  the  negative  entries. 
Moreover,  we  can  assume  that  L  contains  one  negative  and  at  least  two 
positive  diagonal  entries,  for  otherwise  it  would  represent  a  point,  a 
line,  a  plane,  or  a  pair  of  planes,  all  of  whose  intersections  with 
another  quadric  surface  may  be  easily  computed. 

Proceeding  in  this  way,  the  transformed  L  will  have  one  of  the 
following  forms:  Either 

(i)   L   has   three   positive  and  one  negative  entries,  in  which  case  it 

9    o    9    o 

represents  the  sphere  x  +y  +z  -w  =  0;  or 

(ii)  L  has  two  positive  and  two  negative   entries,   in  which   case   it 

9    9    9    9 

represents  the  single-sheeted  hyperboloid  x^+y^-z  -w"  =  0;  or 
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(iii)   L   has   two   positive,  one  negative  and  one  zero  entry,  in  which 

2   2   2 
case  it  represents  the  cylinder  y  +z  -w   =  0  . 


Since  we  have  noted  in  the   preceding   that   all   these  canonical 

surfaces   are   rational,   it   follows   that   every   quadric  surface  Is 
rational.   These  surfaces  have  relatively  simple  structure,  which  makes 

it   easy   to  deal  with  their  intersection  curves.   This  will  be  done  in 

the  present  section,  ultimately  to  yield  a  relatively  simple  technique 

for   the   analysis  of  general  Boolean  combinations  of  bodies  bounded  by 
quadric  surfaces. 

Let  A  and  Z  be  two  quadric  surfaces,  represented  respectively  by 
the  4x4  matrices  L,  S  in  the  manner  described  above.  Applying  the  same 
projective  transformation  to  both  surfaces,  we  may  assume  Chat  L  has 
one  of  the  three  canonical  forms  listed  above,  and  that  S  is  still 
symmetric.  Thus  L  represents  either  a  sphere,  a  single-sheeted 
hyperboloid,  or  a  cylinder,  and  we  are  left  with  the  problem  of  finding 
the  intersection  curve  of  such  a  surface  with  a  second,  arbitrary 
quadric  surface.  This  problem  was  studied  by  J.  Levin  [L] ,  who  gave  a 
fairly  systematic  account  of  these  intersections.  However,  since  the 
situation  is  simpler  when  viewed  projectively  than  when  only  affine 
simplifications  are  used,  we  will  analyze  these  intersections  from  a 
projective  viewpoint.  To  obtain  the  simple  classification  of  the 
intersection  curves  which  we  seek,  we  transform  the  homogeneous  4-space 
projectively  so  as  to  bring  both  quadratic  forms  (Lx,x)  and  (Sx,x)  to 
simple  form,  from  which  the  intersection  curves  can  easily  be  obtained. 

First  we  consider  the  case  in  which  L  is  nonsingular 
(corresponding  to  the  preceding  cases  (i)  and  (ii)),  but  abandon  the 
restriction  that  the  space  X  in  which  L  and  S  act  is  4-dimensional, 
since  it  is  not  hard  to  treat  the  general  n-dimensional  case. 

Our  task  is  to  reduce  the  matrix  S  to  its  simplest  equivalent  form 
while  maintaining  some  equally  simple  form  of  L,  i.e.  to  replace  S  by 
a  simpler  matrix  T'ST  where  T  is  some  real  transformation  for  which 
T'LT  =  L,  or  for  which  T'LT  retains  some  equally  convenient   form.    We 
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will   always   keep   L   in  a   form  satisfying   L^   =   I,   and  then  the 
eigenvalues  of  LS  are  invariant  under  such  transformations  T,  since 

LT'ST  =  (LT)"^ST  =  T"^LST. 

In  what  follows,  we  will  show   that   these   are   essentially   tlie   only 
invariants  of  the  pair  L,  S  . 

Write  [x,y]  for  the  nonsingular  Indefinite  Inner  product  Lx«y. 
Then 

[LSx,y]  =  Sx«y  =  x«  Sy  =  Lx«  LSy  =  [x.LSy], 

so  that  LS  Is  self-adjoint  relative  to  the  bilinear  form  [x,y]. 

The  (generalized)  eigenspace  X^  belonging  to  an  eigenvalue  \  of  LS 
Is  the  set  of  all  vectors  x  such  that  (LS-XI)'^x  =  0  for  some  Integer  n. 
The  space  X  decomposes  Into  the  direct  sum  of  such  elgenspaces,  each  of 
which  is  invariant  under  LS.  It  is  plain  that  for  each  eigenvalue  \i  t 
\,  the  restriction  of  LS  -  XI  to  X^  Is  nonsingular.  Hence  for  each  y  e 
Xy  there  exists  a  z  e  Xj  such  that  (LS  -  XI)'^z  =  y.  Thus  for  each  x  z 
\,  [x,y]  =  [x,(LS-\I)'^z]  =  [(LS-XD^x.z]  =  0.  This  shows  Chat 
elgenspaces  X^^  ,  X^  belonging  to  different  eigenvalues  X  ,^i  are 
orthogonal  to  each  other  In  the  bilinear  form  [x,y]. 

To  reduce  L  and  S  simultaneously  to  the  desired  canonical  form  it 
Is  sufficient  to  analyze  the  restriction  of  LS  to  each  of  the  spaces 
'^\  ,  X  real,  and  to  each  pair  of  spaces  X;^  ,X^,  X  complex.  First 
consider  a  generalized  eigenspace  X^  with  X  real.  In  this  space,  the 
operator  R  =  LS  -  X  I  Is  still  symmetric  relative  to  the  nonsingular 
bilinear  form  [x,y]  and  satisfies  R^=0  for  some  (smallest)  n,  so  that 
R"^"^  *  0.  Thus  there  exist  x  and  y  In  X^  such  that  [y,R"~-^x]  *  0.  This 
Implies  that  there  exists  a  z  e  X-^  such  that  [z,r"~  z]  *  0,  since 
otherwise  x,  y  would  satisfy 

0  =  [(x+y),R""^(x+y)]  =  [y,R""^x]  +  [x,R'^-^y] 
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=  2[y,Rf^-lx]  t    0. 

Suppose  without  loss  of  generality  that  [x,R""^]  f-  0.  Then  the  the 
space  Y^  spanned  by  {x,Rx, . . . ,R"~  x}  is  Invariant  under  R  (I.e.,  under 
LS),  and  the  restriction  of  [x,y]  to  this  space  Is  nonslngular. 
Indeed,  if 

u  =  a  .R Jx  +  a  .^^RJ'^^x  +  .  •  •  +  a^^.^R'^'-^x, 
where  ex  •  *  0  Is  such  that  [u,R"^]  =  0  for  all  m  >  0,  then 

0  =  [u,R""^"Jx]  =  aj[RJx,R""^"Jx]  =  a  ^  [x,r"~^x]  jst  0, 

a  contradiction. 

Since  the  space  Y|^  is  Invariant  under  R  and  since  [x,y]  is 
nonslngular  on  Y^,  Its  orthogonal  complement  Y^'  (relative  to  the  inner 
product  [x,y])  has  these  same  properties.  Thus  the  argument  given  In 
the  preceding  paragraph  can  be  applied  to  Y^^ ' ,  from  which  It  plainly 
follows  that  X^  decomposes  Into  a  direct  sum  of  subspaces  Y^  +  Y9  +  ... 
+  Yj.,  mutually  orthogonal  relative  to  the  Inner  product  [x,y],  In  each 
of  which  the  operator  R  acts  as  a  shift  (i.e.  each  of  these  spaces  has 
a  basis  {x,Rx, . . . ,R"~^x} )  for  which  [x,R^~^x]  *    0. 

Consider  a  single  such  space  Y  =  Y^,  and  let  Its  dimension  be  n, 
so  that  R^y  =  0  for  all  y  e  Y.  By  modifying  x  appropriately  we  can 
assume  without  loss  of  generality  that  [x,Rix]  =  0  except  for  j  =  a-1. 
Indeed,  suppose  that  this  Is  false,  and  consider  the  largest  j<n-l  for 
which  a  =  [x,RJx]  *  0.  Then  for  each  3  we  have 

[x  -  0R"-J-lx,  RJ(x  -  SR'^'J"^)] 
=  [x,RJx]  -B[x,R"-^]  -3[x,Rn-lx] 
=  a  -  2B  [x,R""^]. 

Hence    if   we   put  B    =  a /(2  [x,R"~^] )   and   x'      =     x     -     6R""J"^x,      we     have 
[x',RJx')    =  0.   Moreover,    for  each   j<l<n-l   we   have 
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[x',Rix']  =  [x  -  3R""J"^x,R^(x  -BR""J~^x)]  =  [x.Rix] 

so  that  [x'.Rix']  =  0  for  j<i<n-l  and  [x',Rr>-lx']  ^  0. 

It  is  plain  that  {x' ,Rx' , . . .  ,R'^~^x' }  spans  the  same  space  Y  as 
{x,Rx, .  . .  ,R'^~^x} .  Thus  if  we  replace  x  by  x',  j  is  reduced  by  at  least 
1,  proving  that  there  exists  an  x  as  above,  but  with  [x,RJx]  =  0  for 
all  (K  j<n-l. 

Multiplying  x  by  a  suitable  real  constant,  we  can  suppose   without 

loss   of   generality   that   [x,R'^~^x]   =  c  =  1  or  =  -1.   Note  then  that 

[RJx,R^x]  =  c  if  j+k  =  n-1  but  0  otherwise.    Thus   the   bilinear   form 

[x,y]  has  one  of  at  most  two  canonical  forms  on  the  space  Y  (on  which  R 

has   a    fixed   canonical    form).     Specifically,    in    the   basis 

{x,Rx, . . . ,R"~  x}  the  quadratic  form  [y,y] 

n-1 
(where  we  write  y  =   E    y  R'^  ),  is 

m=0 

F(y)  =  ±(  2yQy^_;^  +  ...  +  Zyj^yn-m-l)*         ^   =  (n-2)/2  if  n  is  even; 
F(y)  =  ±(  '^Yqy^-i    +  ...  +  Zyj^-iy^-xn  +  Ymym^'     '^  =  (n-l)/2  if  n  is  odd, 

Moreover,  in  this  same  basis,  the  form  (Sy,y)  =  [LSy,y]  =  [ (R+a I)y,y ] 
is 

G(y)  =  ±(  2yQyj,_2  +  ...  +  2y^yj,_ni-2)  +  AF(y), 

m  =  (n-3)/2  if  n  is  odd; 

G(y)  =  ±(  2yoy^_2  +  ...  +  2yn,-iyn-m-l  +  Ymym)  +  ^^(y). 

m  =  (n-2)/2  if  n  is  even. 


It  is  clear  that  if  n  is  even  the  signature   of   F   is   (n/2,n/2), 

while   if   n   is  odd  the  signature  of  F  is  either  ((n-l)/2,  (n+l)/2)  or 

((n+l)/2,  (n-l)/2).   The  differing  signatures  distinguish  the  two  cases 

in  which  n  is  odd,  but,  even  if  n  is  even,  the  cases  [xjR^'^-x]  =  1  and 
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[x,R'^~^x]  =  -1  are  easily  seen  to  be  non-isomorphic.  It  Is   convenient 

to  introduce  symbolic  designations  for  these  various  cases;  we  will  use 

the  symbols  [X,m,ra+1]  and  [A ,m+l ,m]  to  designate  the  two  cases  in  which 

n   is  odd,  and  [X,m,m,+],  [X  ,m,m,-]  to  designate  the  two  cases  in  which 
n  is  even. 

Next  consider  a  pair  X,  X"  of  mutually  conjugate  complex 
eigenvalues  of  LS,  and  the  generalized  eigenspaces  X,  ,Xj-  associated 
with  them.  Since  LS  is  a  real  matrix,  it  is  clear  that  x  e  X,-  if  x  e 
X)^  (where  x  denotes  the  complex  conjugate  vector  of  x.  )  The  analysis 
that  we  have  given  applies  without  change  to  the  space  X,  even  if  X  is 
complex  (except  that  for  X  complex  the  vectors  which  appear  in  this 
analysis  will  not  be  real.)  Hence  X-^^  can  be  written  as  a  direct  sum  of 
mutually  orthogonal  subspaces  X^^  =  Y^  +  ...  +  Y.^,  each  of  which  has  a 
basis  {x,Rx,.  .  .  ,R"~^},  where  R  =  LS-Xl,  R"x  =  0,  and  [RJx,R"x]  =  1  if 
j+m=n-l  but  0  otherwise.  Consider  one  of  these  spaces,  e.g.  Yi,  and 
as  a  real  basis  for  the  direct  sum  Y^  +  Y^  take  the  vectors  a^  =  1/7  2 
(Xj  +  x-j)  and  b .:  =  l/i/2  (x^-xj),  where  x.:  =  R-^x,  j=0,...,n-l.  Since 
the  vectors  x^.x^  belong  to  different  eigenspaces,  they  are  orthogonal 
in  the  form  [u,v].   It  follows  easily  that 

[aj,b_^]  =  0,      for  all  j,m, 

'■^j'^m^  =  ~  [t)j,b^]  =  1,     if  j+m  =  n-1 ,  but  0  otherwise. 


rms 


Moreover,   since  LSxj  =  Xj+i  +  Xxj,  and  since  X  can  be  written  in  te 
of  its  real  and  imaginary  parts  as  X  =  a,  +  iao,  we  have 

I^Saj  =  1//2  LS(xj  +  Xj)  =  aj+i  +  1//2  (Xxj  +  Xx.) 
=  a^+i   +  /2  Re(Xxj)  =  aj+i  +  a^aj  -  a2bj  , 

and  similarly 


LSbj  =  bj+i  +  /2  Im(Xxj)  =  bj+i  +  a2aj  +  a  ^b j 
If  the  dimension  of  Yj  is  n,  it  is  plain  that  the  signature  of  [x,y]  on 
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'^-l  =  Yi  +  Y^  is  (n,n).  In  the  basis  {api, .  .  .  ,g.^-i,hQ, . . .  ,h^-i}  for  Zj, 
the  quadratic  form  (Lu,u)  has  the  following  representation  (derived  by 

n-1 
writing  u  =  Z  (^Yj^j  +  Zjbj)  ): 
j=0 

F(y,z)  =  2yQyj^_^  +  ...  +2Y^y^^.i    -ZzQZn-l  "  ...  -2z^z^_^_l , 

m  =  (n-2)/2,   if  n  is  even; 

F(y,z)  =  2yQy^_^  +  ...  +2Y^-iy^-^   +yn,yn,  -2zoZn-i  -  ...  -2z^_iZn-m  "Zm^m. 

m  =  (n-l)/2,   if  n  is  odd. 

Moreover,  in  this  same  basis  (Su,u)  =  [LSu,u]  =  [(R  +  \I)u,u]  has  the 
representation 

G(y,z)    =   2yQy^_2  +   ...    +2yjnyn-m-2   -2zoZn-:    -    ...    -2z^z^-^^2 

+  aiF(y,z)   +  2a2(yoZn-l    +   •••   +  Yn-l^o). 

ra   =    (n-3)/2,      if    n   is    odd; 

G(y,z)    =   2yoyn-2  +   ...    +2yjn_iyn-m-2   ^Ymym  "220^0-2   "    •••    -2zm-l2n-m-2   "Zm^ra 
+  aiF(y,z)   +  2a2(yoZn-l    +   •••    '^  Yn-l^o). 

m   =    (n-2)/2,      if    n   is    even. 

These  pairs  of  quadratic  forms  are  defined  uniquely  by  the  complex 
eigenvalue  X  and  by  the  dimension  n,  and  we  can  therefore  designate 
them  by  [X,n,n],  X  complex. 


The  following  theorem,  which  should  be  compared  to  the 
corresponding  result  for  the  complex  case,  given  in  [HP52,  Vol.  II, 
Chap.   XIII,  Section  10],  summarizes  the  preceding  analysis: 

Theorem  J_:  Let  L  and  S  be  a  pair  of  symmetric  matrices  acting  in  a  real 
vector  space  X,  and  suppose  that  L  is  nonsingular.  Then  X  can  be 
decomposed  as  a  direct  sum  of  canonical  subspaces  Y  which  are 
simultaneously  orthogonal  to  each  other  in  the  two  bilinear  forms 
(Lx,y)  and  (Sx,y).   A  descriptor  characterizing  the  representation  in  a 
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standard  basis  for  Y  of  each  of  these  two  forms  is  associated  with  each 
such  Y.  The  following  table  shows  the  dimension  of  the  space  Y  and  the 
signature  of  the  form  (Ls,y)  on  Y  for  every  possible  descriptor: 

Descriptor  Dimension  of  Y^        Signature  of  (Lx,y )  on  Y 

[X,m,m+1],  X  real  2m+l  (m,m+l) 

[,\,m+l,ra],  X    real  2m+l  (m+l,m) 

]X,m,m,+],  X  real  2in  (m,m) 

[X,ra,m,-],  X  real  Zm  '            (m,m) 

[X ,m,ra] ,  X  complex  2m  (m,m) 

J. 
The  representations  which  the   forms   (Lx,y)   and   (Sx,y)   take   on   in 

standard   bases  for  spaces  having  these  descriptors  are  detailed  in  the 

preceding  paragraphs. 

To  complete  our  analysis  it  remains  to  handle  the  situation  in 
which  L  is  singular.  For  this,  we  have  only  to  carry  over  the  analysis 
described  in  [HP52],  pp. 281-286,  which  we  review  for  the  reader's 
convenience.  As  previously,  let  X  be  the  finite  dimensional  space  on 
which  L  and  S  act.  We  shall  show  that  if  L  is  singular  X  can  be 
decomposed  into  at  least  two  subspaces,  mutually  orthogonal  with 
respect  to  both  L  and  S,  in  one  of  which  the  quadratic  forms  F  =  Lv« v 
and  G  =  Sv  v  have  one  of  the  four  following  representations: 

Case  [2n,2n+l]: 

F(x)  =  2x2X3  +  ...  +  2x2nX2n+l,  n> 0 

G(x)  =  2x^x2  +  ...  +  2x2n-ix2n  ±  ^In+l^ 

(simplest  case:  O.x^) 

Case  [2n-2,2n]: 

F(x)  =  2x2X3  +  ...  +  2x2n-2X2n-l.      n> 1 
G(x)  =  2x^x2  +  ...  +  2x2n-lX2n 
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(simplest  case:  0,2xy) 

Case  [2n-l , 2n] : 

_   F(x)  =  2x2x3  +  ...  +  2x2n-2X2n-l  -  X2n^ »     ^^ 
G(x)  =  2x^x2  +  •••  +  2x2n-ix2n 


(simplest  case:  y-,2xy) 


Case  [2n, 2n] : 

F(x)  =  2x2x3  +  ...  +  2x2nX2n+l.     "> 1 
G(x)  =  2x^x2  +  •••  +  2x2n-lX2n 

(simplest  case:  2xy,2yz) 


We  can  use  the  symbols  [2n,2n+l],  [2n-2,2n],  [2n-l,2n],  [2n,2n] 
respectively  as  designators  for  these  linear  spaces.  The  following 
table  gives  the  dimensions  and  ranks  of  L  and  S. 


Designator 

[2n,2n+l] 
[2n-2,2n] 
[2n-l,2n] 
[2n.2n] 


We  will  use  the  symbol  [0]  to  designate  a  1-dimensional  space  in 
which  both  the  quadratic  forms  L  and  S  are  0, 

Note  that  S  is  nonsingular  in  all  these  spaces  with  the  exception 
of  the  space  [2n,2n].  Thus  each  of  the  spaces  [2n,2n+l],  [2n-2,2n], 
[2n-l,2n]  can  be  regarded  as  a  kind  of  'reverse'  of  spaces  (or  sums  of 
spaces)   that   we  have  already  encountered,  obtained  by  interchanging  L 


Dimension 

Rank 

of 

L 

Rank   of    S 

2n+l,n>0 

2n 

2n+l 

2n,n>  1 

2n-2 

2n 

2n,n>l 

2n-l 

2n 

2n+l,n> 1 

2n 

2n 
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and  S.  By  considering  the  form  of  SL,  it  is  easily  seen  that  in  this 
sense  [2n+l,2n]  is  the  reverse  of  either  [0,n,n+l]  or  [0,n+l,n];  that 
[2n-2,2n]  is  the  sum  of  two  isomorphic  spaces,  each  of  which  is  the 
reverse  of  one  of  the  spaces  [0,k,k+l],  [0,k+l,k],  [0,k,k.+],  or 
[0,k,k,-];  and  that  [2n-l,2n]  is  the  reverse  of  either  [0,n,n,+]  or 
[0,n,n,-].   The  exceptional  space  [2n,2n]  is  its  own  reverse. 

In  view  of  those  facts,  it  is  convenient  to  introduce  the 
notations  [X,m,m+1]  ,...   for  the  reverse  of  [X,m,m+1]... 

To  handle  the  case  in  which  L  can  be  singular,  we  proceed  by 
induction  on  the  dimension  of  X.  If  L  is  nonsingular,  there  is  nothing 
to  prove,  so  that  we  can  assume  that  the  nullspace  N=N|^  of  L,  i.e.  the 
set  \  =  {x  e  X:  Lx=0},  is  nonempty.  If  x« Ly  =  0  for  all  y  then  Lx  =  0 
and  conversely;  hence  the  range  R  =  Rj^  of  L  is  the  orthocomplement  of 
^L,  giving  the  direct  sum  decomposition  X  =  Ny^  +  R^.  If  there  exists  a 
nonzero  x  £  N-j^  such  that  Sx  =  0,  then  we  can  write  X  as  the  sura  of  the 
1-dimensional  space  (x)  and  its  orthocomplement  Y,  and  these  spaces  are 
orthogonal  in  both  forms  L  and  S.  Hence  in  this  case  we  achieve  a 
complete  canonical  decomposition  of  X  by  decomposing  Y  inductively. 
Thus  we  need  only  consider  the  case  in  whch  S  is  nonsingular  on  Nt  .  If 
there  exists  any  x  g  N^  such  that  x- Sx  #0,  we  can  let  Y  be  the 
S-orthocorapleTient  of  the  1-dimensional  space  (x),  and  then  again  X  is 
the  direct  sum  Y  +  (x)  of  two  spaces  orthogonal  in  both  forms  L  and  S, 
so  once  more  we  can  proceed  inductively.  In  all  remaining  cases  we 
will  have  x« Sx  =  0  for  all  x  e  N,  which  implies  that  y« Sx  =  0  for  all 
x,y  e  N,  so  that  S  is  a  1-1  mapping  of  N  into  its  orthocomplement  R. 

Take  some  x^  e  N.  Let  Z  be  the  space  spanned  by  x^  and  Sxq.  The 
form  (u« Sv)  is  nonsingular  in  Z,  since  if  z«S(ax  +  bSx^)  =  0  for  all  z 
e  Z,  then  putting  z  =  x^  we  have  b|Sxo|-  =  0  since  x^. Sxg  =  0,  and 
therefore  b  =  0,  but  now  putting  z  =  Sx^  gives  a  =  0.  It  follows  that  X 
can  be  written  as  a  direct  sum  X  =  Z  +  Z*,  where  Z*  is  the 
S-orLhocomplement  of  Z.  Note  that  Z  and  Z*  are  S-orthogonal ,  but  not 
necessarily  L-orthogonal. 
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Since  x^.  Sxq  =  0,  the  noasingular  form  S  is  neither  positive 
definite  nor  negative  definite  on  the  space  Z,  and  hence  it  must  liave 
signature  (1,1).  Thus  in  appropriate  coordinates  it  has  the 
representation  x--y'-,  and  in  these  coordinates  x  will  appear  as  a 
2-dimensional  vector  (a,b)  satisfying  a' — b^  =  0.  Put  e^^  =  Xq,  and  let 
^2  be  the  vector  with  coordinates  c(a,-b)  in  these  same  coordinates, 
where  c  =  (2a-)"^  then  ej^.  Sei  =  e2*  Se2  =  0  and  ei«  Se2  =  1.  1-et 
e3,...en  be  a  basis  for  Z*.  In  the  basis  for  Z  consisting  of  the 
vectors  e^  and  e2,  the  form  S  has  the  matrix  A2  appearing  just  below. 
Moreover,   the   symmetric  linear  transformation  representing  the  form  L 

satisfies  Le^  =  0.  In  the  following  discussion  we  will  make  use  of  this 
basis,  and  write  u« v  for  a  positive-definite  inner  product  in  which  the 
vectors  e-  are  orthonormal. 

\t  this  point  we  begin  a  subsidiary  induction.  Let  n  be  the 
dimension  of  X,  let  A2  designate  the  2x2  matrix 


[;;] 


and  for  each  even  m  let  A^^  designate  the  mxm  block  diagonal  matrix 
built  from  m/2  conies  of  A2.  For  m  odd,  let  A^  designate  tlie  mxm 
matrix 


^m-1 


Let   B,C   designate  either  L,S  or  S,L.   Suppose  that  for  some  k  we  have 
found  a  basis  in  which  B,C  have  the  respective  forms 


(1) 


^1. 


and 


^k+l 


C 


where  B'  is  (n-k)x (n-k)  and  C  is  (n-k-1 )x (n-k-1 ) .   (Note  that   in   the 
preceding   paragraph,   we  have  in  fact  defined  such  a  basis  for  k  =  1.) 
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Then  we  will  show  that  there  also  exists  a  basis  in  which  C  has   the 
same  form,  but  in  which  B  has  one  of  the  three  following  forms: 


(2) 


1 


or 


\ 


±1 


or 


Ak+2 


where  B"  is  (n-k-l)x  (n-k-l)  and  B'"  is  (n-k-2)x  (n-k-2) .  In  the  first 
two  of  these  three  cases,  there  evidently  exists  a  pair  of  spaces  of 
respective  dimensions  k+1  and  n-k-l  which  are  orthogonal  to  each  other 
in  both  forms  L,S,  and  it  is  also  plain  that,  in  the  first  of  these 
spaces,  the  forms  L,S  have  one  of  the  four  previously  listed  canonical 
representations.  In  the  third  case  we  simply  interchange  R,C,  increase 
k  by  1,  and  continue  our  subsidiary  induction  on  k,  which  must 
eventually  terminate. 

For  the  inductive  step  we  argue  as  follows.  Let  the  basis  in 
which  B,C  have  the  forms  (1)  be  e^  , . . .  ,e,^,ei<.+i , . .  .e^.  '^^  "^'^k+l  =  0  we 
have  the  first  case  (2).  If  not,  modify  the  basis  so  that  e.,..,en, 
where  j  >  k+I ,  is  a  basis  for  the  nullspace  of  ?>' .  This  gives  B  the 
form 


(3) 


Av 


1 


m 


where  0^  designates  an  mxm  square  matrix  of  zeroes,  m  =  n- j+1 ,  and  Bq 
is  nonsingular,  say  of  dimension  (and  thus  rank)  r.  Use  an  auxiliary 
positive  definite  inner  product  in  which  the  basis  ei...e  is 
orthonormal,  and  let  E  and  F  denote  the  orthogonal  projections  of  X 
onto  e^^2.' •  •  .ej_i  and  onto  e^,...,ej_i  respectively.  Then  EBF.  omits 
one  row  and  one  column  of  B^,  and  hence  its  rank  is  either  r-1  or  r-2. 
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Suppose   first  that  the  rank  of  EBE  is  r-1 ,  i.e.   identical  to  the 
rank  of  EB,  so  that  these  two  matrices  have  the  same  range.   Then  there 
exists   a  u  =  Eu  such  that  EB(e^^^^-u)  =  0.  Plainly  ei^+^'u  =  0  =  (Cei^)«u 
since  Ce^   =   e^+[ .   Define  a  linear  transformation  T  on  X  by  putting 

Tx  =  X  -  (e|^+px)u  +  (Cu«x)e{^  -  1 /2(e^+^' x)(Cu«  u)e^. 

Then  if  Tx  =  0,  we  have  e^^^p  Tx  =  e^^+px,  so  Tx  =  x  +  (Cu'x)ei^.  Hence 
0  =  Cu« Tx  =  Cu« X  +  (Cu« x)(Cu« e^)  =  Cu* x,  and  therefore  x  =  0,  showing 
that  T  is  nonsingular.  Plainly  Te^  =  e^^  for  i  =  l...k.  Moreover  FBFej^ 
=  0,  so  that 

(Ty).FBF(Tx)  =  (y-(e^^^>.y)u  +   ((Cu-y)  -  l/ZCe^^+l*  y)(Gu.  u)  )e,J 
•FBF(x-(e^^-^px)u  +  ((Cu-x)  -  1 /2(e,.+i.  x)  (Cu*  u)  )ei^) 
=  ((e^.+l«y)  (e^^+i-u)  +  Ey  )•  FBF(  (e^+1' x)  (e;^+^-u)  +  Ex) 
=yEBEx  +  (e^^py)  (e^+i«  x)  (e^+i    -  u)«  B(e|^+l-u) . 

Thus  in  the  basis  Tej^ , . .  .Te^^,  the  form  B  has  the  representation 
A, 


^o' 


where   B^'  has  one  less  row  and  column  than  B^.   Plainly  c  *  0  since  Bq 
is  nonsingular;  hence  we  can  normalize  and  take  c  =  ±1. 


Since  Ce|^  =  e^^+j^  and  e^^*  Ce^  =  u*  Ce|^  =  e^^*  Cu  =  0,  we  have  also 


Ty.C(Tx)  =  (y-(e(^+py)u  +  ((Cu-y)  -l/2(e,^+l«  yXCu*  u)  )e,^) 
•C(x-(e^^^.x)u  +  ((Cu'x)-l/2(ei^+l'x)(Cu«u 
=  y  Cx-(e,  , ,.  v")u«  Cx  +  (CCu'v)  -1/2  (ei,a.i  •  v)  (Cu 


y  Cx-(e(^^P  y)u«  Cx  +  ((Cu«y)  -1/2  (ej^+p  y)  (Cu«  u)  )e|<.«  Cx 


•l^+l'x)(yCu-(e|^+l«y)u«Cu)+((Cu«x)-l/2(e,^+l«x)(Cu«u))yCe,< 
Cx-(ej^^p  y)u«  Cx  +  (e,^^px)(u'Cy)  -  (ej^^^*  x)  (y  Cu) 
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+  (e^^+|•x)(e^^+l•y)u.  Cu  -l/2(e,^+i- y)e,^^.l«  x)Cu.  u 
=  y  Cx 

In  this  new  basis  B  plainly  has  the  second  of  the  forras  (2),  completing 
our  treatment  of  the  case  in  which  the  rank  of  EBE  is  r-1. 

Next  suppose  that  the  rank  of  EBE  is  r-2,  so  that  there  exists  a 
nonzero  vector  w  e  EX  such  that  EBEw  =  0,  ie.  EBw  =  0.  Change  the 
basis  to  make  w  =  e^^',.  Then  in  (3)  the  submatrix  B^  will  have  the 
form 


(4) 


c  0  .  .  .  0 


^o' 


*  0 


where  the  asterisks  denote  elements  that  may  be  nonzero,  and  B  '  is  an 
(r-2 )x (r-2)  matrix  which  is  plainly  nonsingular.  Define  F  and  E  as 
before.  Since  '^  '  is  nonsingular  there  clearly  exists  a  a  =  Eu 
satisfying  EB(e^^j^-u)  =  ce],^.2.  (The  vector  u  is  determined  by  this 
condition  only  up  to  addition  of  a  multiple  of  e^^^O')  Computing  just  as 
before  we  find  that 

(Ty)' FBF(TX)  =  y  EBEx  + 

c((ek+l-y)  (ek+2-x)  +  (e,^+i-x)  (e^+2-y)) 
+  a(e^+^.y)(e,^4.1.x), 


where  a  =  (e^^.i -u)«  BCe^^ .  i -u) ,  and  we  also  ha 


ive 


Ty.  C(Tx)    =  yCx. 


We      must      have      c     *      0,    since   B^   is    nonsingular;    ch.inging    the   scale    of 
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^v+? .  we  can  therefore  assume  that  c  =  1.  Write  u  =  u^  +  be^,^7 ,  where  b 

is  at  our  disposition,  since,  as  already  noted,  u  is  determined  only  up 

to  addition  of  a  multiple  of  e^^^.^.   Then,  since  Be^^+9  =  ce'^^.^  =  e]._^|  we 
have 

=  ^^k+1  -  ^o)*2(^k+i-^o)   +  -'^(^k+r^o^^k+l 

=  ^^k+1  -  %>-^^^k+l-%>  +  2b. 

Choosing  b   =  ~l/2(eij^^i-u  )»B(ei^^.i-u  )  gives  a  =  0.  Hence  in  the  basis 
Tep...,Tej^,  B  has  the  form 


0  1 

1  0 


where  B''  is  an  (n-k-2)x (n-k-2)  matrix.  This  is  the  third  case  of  (3), 
and  simply  continues  our  induction,  which  must  eventually  terminate 
with  one  of  the  first  two  cases  of  (3). 

This  completes  the  proof  of  the  following  theorem,   which   extends 
Theorem  1  to  the  singular  case: 


Theorem  2:  Let  L  and  S  be  a  pair  of  symmetric  matrices  acting  in  a 
real  vector  space  X.  Then  X  can  be  decomposed  as  a  direct  sum  of 
canonical  subspaces  Y  which  are  simultaneously  orthogonal  to  each  other 
in  the  two  bilinear  forms  (Lx,y)  and  (Sx,y).  A  descriptor 
characterizing  the  representation  in  a  standard  basis  for  Y  of  each  of 
these  two  forms  is  associated  with  each  such  Y.  Each  such  descriptor, 
and   its  associated  space  and  form,  is  either  one  of  those  appearing  in 

■k  "k 

Theorem  1,  or  one  of  the  spaces  described   by   [0,m,m+l]  ,   [0,m+l,m]  , 

"k  k  "k 

[0,m,m,+]  ,  or  [0,m,m,-]  ,  where  D  designates  the  result  of  taking  the 
space  and  pair  of  canonical  forms  described  by  D  and  interchanging  the 
two   forms,   or   is  the  1-dimensional  space  with  both  forms  ^ero,  whose 
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descriptor  Is  [0],  or  is  an  odd-dimensional  space  with  descriptor  [n,n 
and  forms 

x^X2+...+X2n-iX2n  and  '^2'^3"'"' *  •■^2n^2n+l  • 


3.  Application  to  Computation  of  Intersection  Curves 

Now  we  can  apply  the  preceding  analysis  of  the  general 
n-dimensional  case  to  the  4-dimensional  case  required  for  the  geometric 
application  which  interests  us.  We  use  the  notations  introduced  at  the 
beginning  of  the  preceding  section.  For  the  moment  we  assume  that  L  is 
nonsingular.  Suppose  first  that  LS  has  a  nonsimple  real  eigenvalue  X. 
Make  X =0  by  subtracting  XL  from  S,  and  let  n  be  the  dimension  of  one  of 
the  spaces  Y  appearing  in  the  canonical  decomposition  of  the 
generalized  eigenspace  corresponding  to  X ,  so  that  a  equals  either  2, 
3,  or  4.  If  n=2,  then  Y  has  coordinates  in  which  (Lu,u)  and  (Su,u)  have 
the  respective  forms  ± 2zw  and  tz".  Thus  if  we  pass  to  inhomogeneous 
coordinates  for  the  whole  of  our  original  Euclidean  3-space  by  dividing 
by  z,  the  equations  (Lu,u)  =  0  =  (Su,u)  take  on  the  form 

Ql(x,y)  +  2w  =  0  =  Q2(x,y)  +  1 

where  Q^  and  Q2  are  quadratic  forms  in  two  variables.  The  second  of 
these  equations  defines  a  quadratic  plane  curve  to  which  we  can  easily 
give  a  rational  parametrization,  and  then  the  first  equation  gives  the 
w  coordinate  for  the  corresponding  points  on  the  intersection  of  A  and 
Z.  Thus  in  this  case  the  intersection  curve  of  A  and  Z  has  rational 
parametrization. 

Next  suppose  that  n=3.  Then  on  Y  the  forms  (Lu,u)  and  (Su,u)  can 
be  represented  as  ±  (2yw  +  z'^)  and  ±2yz.  Pass  to  inhomogeneous  3-space 
coordinates  by  dividing  by  z,  so  that  (Lu,u)  =  0  =  (3u,u)  takes  on  the 
form 


±x^  +  2yw  +1=0  =  ax^  +  2y 
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In  this  case  the  second  of  these  equations  can  be  used  to  eliminate  y 
from  the  first,  giving  the  equation  of  an  elementary  curve  in  the  x,w 
plane,  and  then  the  second  equation  gives  the  y  coordinate  for  the 
corresponding  points   on   the   intersection   of  A   and  Z  .      Hence  the 

9  9       9 

intersection  curve  of  A  and  E  is  simply  (x,-ax  /2,(l±x  )/(ax'')). 


Finally  suppose  that  n=4.  Then  (Lu,u)  and  (Su,u)  can  be 
represented  as  ± (2xw  +  2yz)  and  ± (2xz  +  y").  Here  we  can  pass  to 
inhoraogeneous  coordinates  by  dividing  by  y,  which  gives  the  equation 


xw  +  z  =  0  =  2xz  +  1 

whose  elementary  parametrization  is  (x,w,z)  =  (1/t,  t  /2,  -t/2). 

It  is  also  possible  that  LS  should  have  complex  eigenvalues, 
either  simple  or  nonsimple.  If  LS  has  a  nonsimple  complex  eigenvalue, 
then  by  subtracting  a  real  multiple  of  L  from  S,  and  multiplying  S  by 
an  appropriate  constant,  we  can  suppose  without  loss  of  generality  that 
this  eigenvalue  is  i.  In  this  case  the  preceding  analysis  shows  that 
the  equations  (Lu,u)  =  0  =  (Su,u)  reduce  to 

9       9 

2(xy  -  zw)  =  0  =  x"  -  z^   +  2xw  +  2yz 

We   can   then  pass  to  inhomogeneous  coordinates  by  dividing  by  z,  which 
gives  the  equations 

xy-w=0=x~-l+  2xw  +  2y 

Using  the  first  equation  to  eliminate  w  from  the  second  and  solving  for 

9      9 
y  yields  y  =  (l-x~)/(2x  +2).   Thus  the  intersection  of  A  and  Z  consists 

of  a  nonsingular  rational  curve. 

Next  suppose  that  LS  has  a  simple  complex  eigenvalue.  If  its 
other  eigenvalues  are  real,  then  either  the  preceding  analysis  or  the 
analysis  given  below  will  apply,  so  the  only  case  that  need  concern  us 
specially   is   that   in  which   LS   has   two  pairs  of  conjugate  complex 
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eigenvalues.  Again,  we  can  suppose  without  loss  of  generality  that  one 
of  these  eigenvalues  is  i,  and  then  the  representation  theorem  derived 
above  implies  that  there  exists  a  pair  of  2-dimensional  spaces 
orthogonal  in  both  forms  (Lu,u)  and  (Su,u),  on  one  of  vv^hich  these  two 
forms  have  the  representations  z  -  w-  and  2zw  respectively.  Thus 
(Lu,u)  =  0  =  (Su,u)  can  be  written  as 

Q^Cx.y)  +  z-  -  w2  =  0  =  Q2(x,y)  +  2zw  , 

where  Qj^  and  O2  are  homogeneous  nonslngular  quadratic  forms  in  the 
variables  x  and  y.  Passing  to  inhomogeneous  coordinates  by  dividing  by 
w,  these  equations  become 

Q^Cx.y)  +  z2  _  1  =  0  =  Q^Cx.y)  +  2z  . 

Eliminating  z  from  the  first  equation  gives  4Q,-fQ2~'^  =  0*  Since  Qi  and 
Q2  describe  two  simple  conjugate  eigenspaces  on  which  LS  has  a  pair   of 

conjugate   complex  eigenvalues   we   can  find  coordinates  in  which  they 

9      o 
appear  as  a(x"  -  y^)  +  2bxy  and  2xy  respectively,  so  that  this  equation 

has  the  form 

0  =  a(x-  -  y-^)  +  2bxy  +  x-y-  -  1  =  (a+y-)x-  +  2bxy  -  (1+ay-) 

If  a  =  0  this  equation  is  solvable  in  rational  terms,  giving  an 
intersection  curve  which  is  a  hyperbola  parallel  to  the  x,y  plane. 
Thus  only  the  case  a  *  0  requires  special  consideration.  In  this  case 
we  can  regard  the  preceding  equation  as  a  quadratic  equation  in  x,  and 
then  its  discriminant  is  A  =  b-y-+(a+y-) ( 1+ay-) .  Suppose  first  that  a 
>  0.  Then  it  is  easy  to  see  that  A  >  0  for  all  values  of  y,  so  that  the 
equation  describes  a  curve  consisting  of  two  infinite  branches,  which 
can  be  parametrized  in  terms  of  y  using  rational  functions  of  y 
together  with  a  single  square  root. 


If  a  <  0,  then  A  is  negative  for  y  =  0  and  for  large  y.  Moreover, 
unless  b  =  0  and  a  =  -1 ,  there  always  exists  a  nontrivial  interval  I 
such  that  A  >  0  whenever  y^  e  I.  Thus  the  curve  defined  by  the  original 
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equatlon  is  a  pair  of  symmetric  ovals  which  can  again  he  parametrized 
using  rational  functions  of  y  and  a  single  square  root. 

In  all  remaining  cases  LS  has  eigenvalues  which  are  all  simple  and 
at  most  one  pair  of  these  eigenvalues  is  complex.  Suppose  first  that 
all  four  eigenvalues  Xj^,...,^^  of  LS  are  real  and  simple.  Then  the 
homogeneous  4-space  decomposes  as  a  direct  sum  Y-^  +  ...  +  Y^  of  four 
corresponding    1-dimensional   orthogonal   eigenspaces.    In   this 

7   9   9   9 

decomposition,   the   forms   (Lu,u)   and   (Su,u)   are  ±x~±y''±z  tw      and 

9      9      9      2 

iX  ^x"^X  2y~tA -jz~±X  ^w  respectively.  Subtracting  an  appropriate  multiple 
of  L  from  S  and  passing  to  inhomogeneous  coordinates  by  dividing  by  w, 
we  can  bring  S  to  the  form 

9     2 

ay~  +  bz   =  c 

so  that  the  intersection  of  A  and  Z  is  the  intersection  of  A  (which  is 
either  a  sphere  or  a  hyperboloid)  with  an  elliptic  or  hyperbolic 
cylinder  parallel  to  the  x-axis. 

Finally,  suppose  that  LS  has  two  real  simple  eigenvalues  and  one 
complex  conjugate  pair  of  eigenvalues.  We  can  assume  without  loss  of 
generality  that  the  complex  eigenvalues  are  ±i.  Hence  we  can  decompose 
the  4-space  into  the  direct  sum  Z  +  Yj^  +  Y2 ,  where  Y|,Y2  are  the 
1-dimensional  eigenspaces  corresponding  to  the  two  real  eigenvalues 
X^,X2j  and  where  Z  is  a   two-dimensional   space   corresponding   to   the 

eigenvalues  ± i.   In  this  decomposition  the  forms  (Lu,u)  and  (Su,u)  are 

9999  99 

x^-y-t  z-±  w-  and  2xy±X  iz^X  ^w"  respectively.   Subtracting  an  appropriate 

multiple   of   L   from  S   and   passing   to  inhomogeneous  coordinates  by 

dividing  by  w,  we  can  bring  S  to  the  form 


9      9 

a(x  -  y)  +  2bxy  =  c  , 


so  that  the  intersection  curve  of  A  and  E  is  the  intersection  of  A  with 
a  hyperbolic  cylinder  parallel  to  the  z-axis. 
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Tn  each  of  the  last  two  cases,  the  coordinate  system  that  we  liave 
introduced  represents  the  intersection  curve  as  the  intersection  of  the 
unit  sphere  or  the  single-sheeted  hyperbolold  x  +y~-z'"  =  l  with  a 
vertical  elliptic  or  hyperbolic  cylinder  (or  with  one  or  two  planes) 
whose  base  quadratic  curve  Q  lies  in  the  x,y  plane  and  has  center  at 
the  origin.  First  consider  the  case  in  which  A  is  a  sphere.  Depending 
on  the  relationship  of  Q  to  the  unit  circle,  this  intersection  curve 
can  either  consist  of  one  or  two  disjoint  closed  loops,  four  branches 
meeting  at  two  symmetrically  situated  critical  points,  or  two  symmetric 
isolated  points.  All  these  curves  can  be  conveniently  (though  not 
quite  rationally)  parametrized  by  introducing  a  rational 
parametrization  x  =  x(t),  y  =  y(t)  for  the  curve  Q,  and  noting  that  z  = 
(l-x2(t)-y2(t))l/2^  ,■   ■  • 

Next  consider  the  case  in  which  A  is  the  single-sheeted 
hyperbolold  x''+y"-z-=l.  If  the  cylinder  base  Q  is  an  ellipse,  the 
intersection  can  be  either  a  circle,  a  symmetrical  pair  of  closed 
loops,  four  branches  meeting  at  two  symmetrically  situated  critical 
points,  or  two  symraetric  isolated  points.  If  Q  is  a  hyperbola,  the 
intersection  will  be  either  a  set  of  four  symetrically  situated 
infinite  curves  symmetric  in  the  x,y  plane,  or  two  sets  of  four  such 
curves  meeting  at  two  diametrically  opposite  critical  points.  These 
curves  can  be  parametrized  by  introducing  a  rational  parametrization  x 
=  x(t),  y  =  y(t)  for  Q,  and  noting  that  z  =  (x-(t  )+y-(t  )-l )  ^'''^. 

With  the  exception  of  the  spaces  with  designators  [0]  and  [2,2], 
the  quadratic  forms  corresponding  to  aS  +  6L  are  nonsingular  for  all 
but  finitely  many  pairs  (a  ,3  )  in  every  canonical  space  appearing  in  the 
preceding  discussion.  But  when  there  exists  (a  ,B  )  such  that  aS  +  gL  is 
nonsingular,  we  can  analyze  the  intersection  of  E  and  A  in  the  manner 
explained  in  the  preceding  paragraphs  by  forming  the  intersection  of  A 
with  the  nonsingular  surface  defined  by  (aS-(-6L)u«u  =  0.  Thus  only  cases 
in  which  the  designators  [0]  or  [2,2]  appear  require  special 
consideration.  In  the  first  of  these  cases,  Z  and  A  can  be  represented 
simultaneously  as  vertical  cylinders  with  quadratic  curves  as  bases, 
and  their  intersections  are  simply  vertical  lines.   On  the  other  hand. 
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if   the  space  with  designator  [2,2]  appears,  then  the  Intersection  of  A 
and  Z    can  be  represented  as  the  set  of  points  satisfying 

±x^  +  2yw  =  0  =  ax^  +  2zw. 

Dividing  by  w  to  pass  to  inhomogeneous  coordinates  we  find  see  that  the 
Intersection  is  just  the  rational  plane  curve  (x,±x  /2,-ax  /2). 


4.  Boolean  Combinations  of  Quadric  Surfaces. 

For  a  solid  geometry  modeling  system  one  will  need  to  use  the 
detailed  information  concerning  quadric  surface  intersections  developed 
in  the  preceding  section  to  represent  the  geometry  of  an  arbitrary 
boolean  combination  B  of  volumes  bounded  by  such  surfaces.  The  solid  B 
is  formed  by  applying  intersection,  symmetric  difference  and  union 
operators  to  halfspaces  H^,  each  defined  by  an  equation  of  the  form 

Pi(x,y,z)  >    0  or   ?^(x,yyz)  <    0 

where  each  P-  is  a  quadratic  polynomial.  Let  Z^  be  the  quadric 
boundary  surface  defined  by  P^(x,y,z)  =  0.  If  necessary,  we  can  apply 
infinitesimal  translations  to  the  surfaces  E-  to  ensure  that  they  are 
in  general  position,  i.  e.  that  any  two  such  surfaces  which  meet 
intersect  transversally.  Thus  the  boundary  of  the  solid  B  is  the  union 
of  closed  subregions  R^  of  Z  ^,  where  R^  is  the  intersection  of  E  j^  with 
the  boundary  of  B.  In  a  geometric  modeling  system  which  admits  general 
quadric  surfaces  to  its  vocabulary  of  basic  shapes,  we  will  have  to 
determine  the  precise  geometry  of  the  region  R^  for  each  quadric 
surface  Z  ^. 

To  do  this,  we  consider  a  fixed  surface  Z  and  study  tlie 
intersection  curves  C^  of  Z  and  E^  (i^l).  This  network  of  curves 
divides  E  into  a  collection  of  2-dimensional  cells  k..,  each  of  which  is 
bounded  by  finitely  many  piecewise  dif f erentiable  arcs  ct  lying  on 
various   curves  C^,   the  arcs   lying  along  C^^  being  separated  by  the 
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intersecClon  points  of  C^  with  other  curves  C::.  The  region  R  =  R^  we 
seek  is  a  union  of  certain  of  the  cells  k...  We  shall  show  how  to 
determine  the  geometry  of  the  region  R  by  dividing  it  into  patches, 
each  a  union  of  one  or  more  cells.  Intuitively  speaking,  these  patches 
are  the  connected  subregions  of  R.  More  precisely,  patches  which  meet 
only  at  isolated  points  will  be  considered  to  be  distinct,  i.e. 
'patches'  are  defined  as  closures  of  the  connected  components  of  the 
interior  of  R.  Thus  distinct  patches  have  disjoint  interiors  and  either 
are  disjoint  or  intersect  at  finitely  many  boundary  points. 

We  now  show  how  algorithms  for  manipulating  algebraic  numbers, 
like  those  reviewed  in  [SS82],  can  be  used  to  determine  the  ordering 
along  a  fixed  parametrized  intersection  curve  C.  of  this  curve's 
intersections  with  all  other  curves  C-.  (To  avoid  discussion  of 
'exceptional'  coincidences,  it  is  assumed  in  the  following  discussion 
that  where  necessary  we  have  applied  an  infinitesimal  perturbation  to 
these  curves  to  put  them  into  general  position,  so  that  any  curves 
which  intersect  meet  transversally. )  Recall  from  the  previous  section 
that  each  curve  C^  can  be  parametrized  by  coordinate  functions  Xj(u), 
yj_(u),  and  z^(u),  each  of  which  is  a  rational  combination  of 
polynomials,  or  square  roots  of  polynomials,  in  the  real  parameter  u. 
This  remains  true  even  after  we  apply  projective  transformations  to  the 
parametrizations  given  in  the  last  section  in  order  to  transform  all  of 
them  into  a  common  coordinate  system.  Let  S  and  S.  be  the  quadratic 
forms  which  represent  the  surfaces  I  and  Z^  respectively.  Examination 
of  the  intersection  parametrizations  listed  above  shows  that  except  for 
the  case  in  which  the  product  SS^  has  two  pairs  of  complex  eigenvalues, 
each  such  coordinate  function  either  has  the  form 

p(u)  +  d(u)^/^ 
1  +  u^ 

where  p  and  d  are  of  degrees  at  most  two  and  four  respectively,  or  a 
similar  form  with  denominator  u^,  or  is  a  rational  fraction  with 
denominator  u  and  numerator  at  most  quadratic,  or  is  a  rational 
fraction  with  denominator  l+u^  and  numerator  at  most  biquadratic,  or  is 
a   simple  quadratic  polynomial.   In  the  above  formula,  the  polynomial  d 
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which  appears  under  Che  radical  is  the  same,  up  to  a  constant  factor, 
for  all  three  coordinate  functions  x,  y,  and  z.  The  intersection  points 
of  C.  and  some  other  curve  C.  are  the  same  as  the  Intersections  of  C. 
and  the  surface  Z  ^,  and  may  be  determined  by  substituting  the 
coordinate  functions  above  into  the  defining  polynomial  P.(x,y,z)  =  0 
of  the  surface  E  •.  This  yields  a  polynomial  equation  q^(u)  =  0,  of 
degree  at  most  eight. 

In  the  exceptional  case  that  the  product  SS.  has  two  pairs  of 
complex  eigenvalues,  one  of  the  coordinate  functions  is  of  the  form 
shown  above,  but  with  numerator  multiplied  by  u.  In  this  case,  which 
can  only  occur  when  both  quadratic  forms  have  signature  (2,2), 
substitution  of  the  intersection  parametrizations  in  the  polynomial  P. 
may  yield  a  polynomial  of  degree  twelve.  If  many  such  cases  occur  in  a 
particular  application,  it  may  be  desirable  to  use  the  following 
somewhat  different  approach  so  as  not  to  have  to  work  with  polynomials 
of  degree  more  than  eight.  (This  approach  may  in  fact  be  superior 
whenever  a  hyperboloid  is  dealt  with.)  Choose  a  basis  in  which  the 
surface  E  is  given  by  the  formula  x  =  yz  in  inhomogeneous  coordinates. 
Then  the  intersection  of  E  and  any  other  surface  E-  can  be  obtained  by 
substituting  yz  for  x  in  the  equation  of  E-,  yielding  an  equation  of 
the  form 

fiCy.z)  =  Pi(z)y-  +  Qi(z)y  +  RiCO  =  o 

where  the  coefficient  polynomials  are  at  most  quadratic  in  z.  This 
makes  it  clear  that  the  intersection  curve  can  be  parametrized  by  z,  at 
most  one  square  root  being  involved.  Two  such  curves  intersect  if  and 
only  if,  for  some  real  z,  the  polynomials  f^  and  f  •,  viewed  as 
polynomials  in  y,  have  a  common  real  root.  They  have  a  common 
(possibly  complex)  root  if  and  only  if  their  resultant,  which  is  the 
determinant  of  the  4x4  matrix 
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of  quadratic  polynomials  in  z,  vanishes.  Mote  that  this  resultant  is 
of  degree  at  most  eight.  One  can  eliminate  nonreal  roots  y  by 
determining  the  sign  of  the  discriminant  of  f.,  which  is  only  of  fourth 
degree,  at  each  root  of  the  resultant. 

It  follows  that  the  intersection  points  of  the  fixed  parametrized 
curve  C.  and  the  various  other  curves  C-  have  coordinates 
(x^  (t )  ,y^(t )  ,z  j^(t ) ) ,  where  t  is  a  real  root  of  some  qy  In  order  to 
determine  the  sequence  of  arcs  into  which  the  collection  of  all  such 
intersection  points  divides  the  curve  C.,  we  can  make  use  of  algorithms 
like  those  described  in  the  Appendix  of  [Sb82],  which  perform  exact 
computations  with  algebraic  numbers,  i.e.  roots  of  arbitrary  rational 
polynomials.  Given  a  polynomial  q,  these  algorithms  show  how  to  find 
disjoint  intervals  on  the  real  line,  each  having  rational  endpoints  and 
containing  precisely  one  root  of  q.  By  combining  such  root-isolating 
intervals  corresponding  to  different  polynomials  q-,  one  may  determine 
the  relative  ordering  of  the  collection  of  roots  of  all  the  polynomials 
q.  of  our  discussion.  Thus  we  may  assume  that  each  intersection  curve 
C-  has  been  divided  into  arcs  a,  each  bounded  by  intersection  points  of 
C^  with  other  curves,  and  that  the  ordering  of  these  arcs  and  points  Is 


Our  next  step  is  to  determine  how  the  various  cells  into  which  the 
network  of  curves  C^^  divides  the  surface  Z  combine  to  form  surface 
patches  on  Z  .  To  do  this,  we  map  Y.  to  the  u,v  plane  by  using  the 
inverse  of  the  parametrization  of  the  surface  E,  which  corresponds  to 
one  of  the  three  canonical  diagonal  quadratic  forms  enumerated  at  the 
beginning  of  Section  2  above.  For  instance,  if  Z  is  the  unit  sphere, 
the  inverse  of  the  parametrization 

X  = '-^—^  ;    y  = ^-1-^  ;   z  =  ^  "  -'  "  ^' 


2?  oo'  90 

1  +  u   +  v"  1  +  u*-  +  v^  1  +  u   +  V- 

is  given  by 

u  =  x/(l+z)  ;  V  =  y/(l+z) 

This   transforms   the   network  of   curves   on  E   to  a   topologically 


-29- 

equivalent  network  in  the  plane,  whose  natural  orientation  facilitates 
determination  of  the  structure  of  the  patches  P..  From  now  on,  we 
shall  therefore  treat  the  collection  of  curves,  arcs,  intersection 
points,  and  patches  as  lying  in  the  u,v  plane  rather  than  on  the 
surface  E . 

For  each  arc  a ,  it  is  first  necessary  to  determine  whether  the 
cell  immediately  to  the  left  (or  right)  of  this  arc  lies  in  the  region 
R,  i.e.  on  the  boundary  of  the  solid  B.  To  do  this,  find  a  rational 
number  s  such  that  the  point  p=  (u(s),v(s))  lies  on  a.  Working  in  the 
u,v  plane  and  using  polynomial  calculations  only,  we  can  determine  the 
equation  of  the  normal  line  to  the  arc  a  through  the  point  p.  Nov; 
choose  a  point  p'  inf initesimally  to  the  left  (or  right)  of  p  on  this 
normal  line.  The  cell  containing  p'  is  a  subset  of  the  solid  B  if  and 
only  if  an  appropriate  Boolean  combination  of  conditions 
P^(x(p' ) ,y(p' ) ,z(p ' ) )> 0  is  satisfied,  where  the  P^  are  the  defining 
polynomials  of  the  bounding  surfaces  E.  of  the  solid  B.  In  order  to 
exclude  cells  which  are  interior  to  B,  we  perform  a  similar  test  on  two 
points  q  which  lie  close  to  E  on  the  normal  line  to  E  passing  through 
o'.  If  both  points  tested  lie  in  B,  the  cell  containing  p  is  interior 
to  the  solid  B,  hence  does  not  lie  on  R.  Otherwise,  precisely  one  of 
the  points  q  lies  in  B,  hence  the  cell  lies  in  R  and  on  the  boundary  of 
3.  Note  that  at  least  one  point  q  must  lie  in  B,  for  otherwise  B  would 
contain  a  degenerate  two-dimensional  component;  this  possibility  is 
ruled  out  by  our  prior  assumption  that  the  various  surfaces  E. 
intersect  only  transversally. 

Thus  we  may  assume  that  the  curves  C.  in  the  u,v  plane  are  divided 
into  sequences  of  arcs,  and  that  the  region  to  the  left  of  each  of 
these  arcs  or  to  its  right  (but  not  both)  is  specified  as  lying  on  the 
boundary  of  the  solid  B.  From  this  information,  we  wish  to  reconstruct 
the  global  geometry  of  the  region  R  in  the  u,v  plane  which  represents 
the  intersection  of  S  with  the  surface  of  B. 
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In  order  to  do  this,  divide  the  collection  of  arcs  a  into  a  set  of 
closed  loops  1-  which  bound  the  patches  of  R.  The  patches  so  obtained 
will  not  be  disjoint,  but  nondis joint  patches  will  meet  only  at  a 
finite  collection  of  vertices  which  lie  at  the  intersection  of 
different  curves.  Each  simply-connected  patch  will  be  surrounded  by  a 
single  boundary  loop,  but  patches  with  n  holes  will  have  n+1  boundary 
loops.  Thus,  to  understand  the  topological  situation,  we  need  to  know 
whether  a  given  patch  lies  to  the  left  or  the  right  of  each  of  its 
(oriented)  boundary  loops,  and  must  also  determine  the  manner  in  which 
different  loops  are  nested. 

For  this,  we  can  proceed  as  follows.  Let  A  denote  the  collection 
of  arcs  a..  Begin  by  discarding  from  A  all  arcs  which  are  interior  to 
R.  These  are  detected  by  the  fact  that  points  near  the  arc  on  both 
halves  of  its  normal  line  lie  in  R.  Next,  recall  that  each  arc  is 
oriented  by  the  parametrization  of  the  curve  on  which  it  lies,  so  that 
we  may  legitimately  refer  to  the  tail  and  head  of  the  arc.  (Of  course, 
we  will  have  to  reorient  some  of  the  arcs  in  order  to  link  them  into 
the  loops  which  we  seek  to  construct.) 

Begin  the  algorithm  by  picking  an  arbitrary  arc  cti.  Assume  that 
the  region  to  the  right  of  this  arc  lies  on  R;  if  not,  reverse  the 
orientation  of  the  curve  on  which  the  arc  a,  lies.  The  first  loop  will 
begin  with  a  p  we  will  find  successive  arcs  which  border  the  patch  to 
the  right  of  a  j^  and  add  these  arcs  to  the  loop  which  we  are 
construct  ing. 

For  this,  choose  the  successor  a2  of  the  initial  arc  a^  from  among 
the  collection  of  arcs  in  A  whose  tail  or  head  is  the  head  h,  of  ai 
Reorient  a  2  so  that  its  tail  is  hp  Since  R  lies  to  the  right  of  a^ 
a2  will  be  that  arc  which  is  first  encountered  when  we  rotate  a, 
counterclockwise  while  holding  its  head  fixed.  Note  that  the  clockwise 
ordering  of  the  arcs  with  endpoint  h-^  can  be  determined  precisely  from 
the  parametrizations  (u^(t ) , v^ (t ))  of  the  intersection  curves  C^:  the 
slope  of  each  curve  at  h^  is  the  quotient  of  derivatives  of  u  and  v 
with   respect   to  t,   evaluated  at  a  suitable  algebraic  number  t,  and 
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hence  these  slopes  can  be  ordered  by  the  polynomial  algorithms 
described  in  [  SS82J  .  Once  the  ordering  of  the  slopes  is  known,  it  is 
easy  to  determine  the  clockwise  order  of  the  arcs  around  h,.  No  two 
slopes  can  be  equal,  since  all  intersections  have  been  made 
transversal. 

This  procedure  for  finding  the  successor  to  an  arc  on  a  boundary 
loop  of  R  should  be  iterated  until  the  initial  arc  ai  is  again 
encountered.  At  this  point  we  can  delete  the  collection  of  arcs  on 
this  loop  from  the  set  of  arcs  A.  Then  we  pick  a  new  initial  arc  from 
A,  and  continue  in  this  manner  until  A  is  exhausted  and  the  complete 
collection  of  loops  which  bound  patches  of  R  has  been  constructed. 
Note  that  during  its  construction,  each  loop  has  been  oriented  so  that 
the  patch  of  R  which  it  bounds  lies  its  right.  It  remains  to  determine 
whether  the  patch  lies  inside  or  outside  this  loop;  for  this  we  need  to 
establish  the  sequence  of  nesting  of  the  various  loops. 

In  order  to  do  this,  consider  a  collection  of  loops  1-  in  the  u,v 
plane  and  define  the  distance  from  infinity  d.  of  each  loop  1-  as 
follows.  If  1^  is  an  outermost  loop,  define  d.  to  be  1.  Inductively, 
if  1^  is  contained  inside  a  loop  1-  whose  distance  from  infinity  is  n, 
but  if  no  loop  lies  between  1^  and  1-,  define  d.  to  be  n  +  1.  Note  that 
1^  runs  clockwise  if  and  only  if  d.  is  odd  and  that  each  patch  in  R  is 
bounded  by  one  clockwise  (outer)  boundary  loop  and  by  a  number  of 
clockwise  (inner)  boundary  loops  equal  to  the  number  of  holes  in  the 
patch. 

In  order  to  compute  the  distances  d-  for  the  collection  of  loops 
1^,  proceed  as  follows.  Choose  a  point  p  on  a  loop  1  and  let  r  be  an 
infinite  ray  with  endpoint  p.  Assume  that  r  is  in  general  position  with 
respect  to  the  loops  l^,  i.e.  that  at  any  intersection,  r  crosses 
non-tangentially  between  the  outside  and  inside  of  1^^.  Of  course,  the 
number  of  intersections  of  this  ray  with  loops  1,  is  finite,  although 
the  ray  may  intersect  a  given  loop  several  times.  Order  the  list  of 
intersections  of  r  with  the  loops  so  that  the  first  intersection  point 
lies   on  an   outermost  loop.   Define  integers  m.  by  the  condition  that 
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1  is  the  loop  which  contains  the  i-th  intersection  point  on  the  ray. 
Thus  d  is  1.  Let  L  be  the  ordered  list  of  integers  m^ .  Construct  a 
sequence  of  sublists  P.  of  L,  beginning  with  ?^  =  [mj^],  which  will 
record  the  passage  of  a  point  q  traversing  the  loop  diagram  along  r,  as 
follows.  Given  P.,  let  P^+j^  be  the  list  formed  from  P^  by  deleting  the 
last  element  m^  from  P^  if  m^^i^  equals  m;:;  otherwise  by  appending  m.^j^ 
to  P..  In  the  second  case,  record  the  fact  that  Km.^^)  is  nested 
immediately  inside  l(m-),  and  set  dCm^^j^)  to  d(m.)  +  1  ,  if  it  has  not 
yet  been  defined.  The  first  case  occurs  when  q  exits  l(ra.),  the  second 
when  q  moves  deeper  into  the  nest  of  loops.  This  procedure  concludes 
by  setting  the  last  list  Pi„„^  to  the  empty  list  'when  q  exits  the 
outermost  loop.  ..  r-j   ',.. 

If  loops  remain  whose  distance  from  infinity  has  not  yet  been 
calculated,  draw  a  ray  from  one  of  these  loops  and  repeat  the 
algorithm,  continuing  until  d(l.)  has  been  calculated  for  all  loops. 
These  distances,  together  with  the  nesting  information  which  was 
recorded  during  the  algorithm,  completely  determine  the  geometry  of  the 
region  R,  and  hence  determine  the  intersection  of  the  boundary  of  the 
solid  B  with  the  quadric  surface  E.  This  geometric  information  can  be 
used  for  various  computations  involving  the  boundary  of  R,  ^.g.  in 
connection  with  a  drawing  or  surface  area  algorithm.  A  geometric 
modeling  system  supporting  these  applications  is  currently  being 
programmed. 
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